Usenet 1993 July/InfoMagic USENET CD-ROM July 1993.ISO/sources/misc/volume11/ephem4.12/part01

Text File | 1990-03-12 | 48.9 KB | 1,701 lines

Newsgroups: comp.sources.misc
subject: v11i027: repost of ephem, v4.12, file 1
From: [e-mail]@ncs-med.UUCP (Elwood C. Downey)
Sender: [e-mail] (Brandon S. Allbery - comp.sources.misc)

Posting-number: Volume 11, Issue 27
Submitted-by: [e-mail]@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem4.12/part01-repost

[I tested it here on uunet; it unshars fine.

Between this and the last posting (and the botched unofficial patch posting),
I just grandfathered the old single-post style in my poster.  How do software
houses with *production* programs cope?  ;-)  ++bsa]

Several sites are reporting that the first file in my 7-part posting of
ephem was munged. they all reported the same munging, so evidently
it got posted badly. Here it is again; could you please repost it asap?

--------------------------------------------------------------------------
# This is the first line of a "shell archive" file.
# This means it contains several files that can be extracted into
# the current directory when run with the sh shell, as follows:
#    sh < this_file_name
# This is file 1.
echo x Manifest
sed -e 's/^X//' << 'EOFxEOF' > Manifest
XFiles shipped as the ephem v4.12 package:
X
XMakefile    just say "make"
XMan.txt        manual, ready for a very generic printer.
XManifest    this file.
XReadme        check here for hints before building.
Xephem.cfg    sample configuration file. you'll want to edit this for your
X        particular location and interests after reading the manual.
X
XThe source files:
Xaa_hadec.c    convert between alt/az and hour angle/dec.
Xaltmenus.c    draws the three alternate lower screens.
Xanomaly.c    compute anomaly.
Xastro.h        unit conversion macros and planet defines.
Xcal_mjd.c    converters to and from modified julian date.
Xcircum.c    main "astronomy" entry point that finds where anything is.
Xcircum.h    defines Now and Sky structures.
Xcomet.c        compute comet position from elements.
Xcompiler.c    compile and execute general expressions with screen fields.
Xeq_ecl.c    convert between equitorial and eclipitic coords.
Xflog.c        logs and retrieves screen locations for logging purposes.
Xformats.c    basic date, time, prompts, etc formats.
Xio.c        basic cursor and screen io control.
Xmain.c        main loop.
Xmainmenu.c    draws the top screen.
Xmoon.c        compute moon position.
Xmoonnf.c    compute new and full moon dates.
Xnutation.c    compute nutation correction.
Xobjx.c        set and make use of all the object-x info.
Xobliq.c        compute obliquity.
Xparallax.c    functions to compute earth rim parallax correction.
Xpelement.c    basic planet position polynomial coefficients.
Xplans.c        use polynomials to find planet location at any certain time.
Xplot.c        set fields for and manage plots.
Xpopup.c        handle the one-liner "popup" menus at the top of the screen.
Xprecess.c    compute precession correction.
Xreduce.c    convert elliptical elements from one epoch to another.
Xrefract.c    atmospheric refraction model.
Xriset.c        find basic rise/set sideral times of a fixed object.
Xriset_c.c    iteratively solve for local rise/set times of moving objects.
Xscreen.h    define all field locations and codes, extra planet codes.
Xsel_fld.c    handle cursor movement commands in the various screens.
Xsex_dec.c    convert between sexagesimal and decimal notation.
Xsrch.c        set and manage the various search functions.
Xsun.c        compute location of sun at any time.
Xtime.c        manage setting and getting the time from the os.
Xutc_gst.c    convert between UT1 and Greenwich sidereal time.
Xversion.c    current version notice, and revision history comments.
Xwatch.c        manage the screen during sky and solar system displays.
EOFxEOF
len=`wc -c < Manifest`
if expr $len != 2381 > /dev/null
then echo Length of Manifest is $len but it should be 2381.
fi
echo x Readme
sed -e 's/^X//' << 'EOFxEOF' > Readme
XGeneral Notes for version 4.12:
X
X1) Read the generic-printer-ready manual in Man.txt.
X
X2) Ephem does not gracefully handle the fact that there is no year 0 in the
X   Christian calendar. It is best to avoid entering a year of 0 for now.
X
XChange Highlights:
X
X1) changes from 4.10 to 4.12:
X    major change is support for optional object to be specified as either a
X    fixed object, or via heliocentric elliptical or parabolic orbital
X    elements. a simple logarithmic magnitude model is also supported.
X    fix slight bug related to nutation.
X    plot fields in the same order they were selected.
X    add copywrite notice in main.c.
X    note <sys/time.h> in time.c for BSD 4.3.
X    add confirmation question before quitting.
X    add d,o,z special speed move chars.
X    watch: add LST to night sky and maintain RTC time back in main.
X    fix bug in timezone related to daytime flag.
X    add w (week) watch advance key code.
X    allow for precession moving dec outside range -90..90.
X    fix bug that wiggled cursor during plotting in rise/set menu.
X    fix bug preventing DAWN/DUSK/LON from being used in search function.
X    add PC_GRAPHICS option in mainmenu.c (no effect on unix version).
X    fix twilight error when sun never gets as low as -18 degs.
X    always find alt/az from eod ra/dec, not from precessed values.
X    lastmjd in plans.c was an int: prevented needless recalcs.
X
X 2) changes from 4.8 to 4.10:
X    fix transit times of circumpolar objects that don't rise.
X    fix plotting search function when searching is not on.
X    fix Objx rise/set bug.
X    don't erase last watch positions until computed all new ones.
X    added USE_BIOSCALLS to io.c: Doug McDonald's BIOS calls
X    allow dates to be entered as decimal years (for help with plots).
X    remove literal ESC chars in strings in io.c.
X    provide two forms of non-blocking reads for unix in io.c.
X    take out superfluous ESC testing in chk_arrow().
X    guard better against bogus chars in sel_fld().
X    use %lf in scanf's.
X    command-line arg PROPTS+ adds to settings from config file.
X    change (int) casts in moduloes to (long) for 16bit int systems.
X
XIf you have source, here is how you can build ephem:
X
X1) io.c:
X   define either UNIX or TURBO_C in io.c depending on whether you wish to
X   compile under unix or dos/turbo/lattice/microsoft.
X
X   Also in io.c and if you use unix, you have two choices of methods for doing
X   non-blocking reads depending on whether you define USE_NDELAY; use whichever
X   works for you.
X
X   MSDOS users can do cursor control with direct BIOS calls (the default),
X   or with ANSI.SYS by defining USE_ANSISYS in io.c.
X
X2) time.c:
X   Select from two methods of dealing with time from the operating system
X   with the TZA/TZB defines in time.c. If you get link undefines related to
X   time functions try changing to the other form.
X
X   On Sun systems running OS 4.0.3 (or BSD 4.3) or Apollo SR 10.1 use TZB and
X   change <time.h> to <sys/time.h>.
X
X3) mainmenu.c:
X   if you are compiling for an IBM-PC then try #define PC_GRAPHICS for a bit
X   nicer looking way to draw the screen boundry lines.
X
X4) beware that I have not used string.h or strings.h. if your library's
X   strlen() and str.*cmp() functions don't return int (such as long), then you
X   will have to hand add string.h or your own extern declarations. strcmp() is
X   used in io.c; strncmp() in compiler.c and main.c; and strlen() is used in
X   compiler.c, io.c, main.c, mainmenu.c, popup.c and version.c.
X
X5) Ephem can now be built by simply compiling all the .c files and linking them
X   all together. On Unix systems, you must also link with the termcap library
X   (-ltermcap) and possibly the auxiliary math library (-lm) if your default C
X   library does not include all the required transcendental functions.
X
XThe following files are pretty much just pure transliterations from BASIC
Xinto C from machine-readable copies of the programs in Duffett-Smith's book.
XThey have nothing to do with the rest of ephem so they may be used for
Xcompletely different applications if so desired.
X   aa_hadec.c anomaly.c astro.h cal_mjd.c comet.c eq_ecl.c moon.c moonnf.c
X   nutation.c obliq.c parallax.c pelement.c plans.c precess.c reduce.c
X   refract.c riset.c sex_dec.c sun.c utc_gst.c
EOFxEOF
len=`wc -c < Readme`
if expr $len != 4245 > /dev/null
then echo Length of Readme is $len but it should be 4245.
fi
echo x Makefile
sed -e 's/^X//' << 'EOFxEOF' > Makefile
XEPHEM=  aa_hadec.o \
X    altmenus.o \
X    anomaly.o \
X    cal_mjd.o \
X    circum.o \
X    comet.o \
X    compiler.o \
X    eq_ecl.o \
X    flog.o \
X    formats.o \
X    io.o \
X    main.o \
X    mainmenu.o \
X    moon.o \
X    moonnf.o \
X    nutation.o \
X    objx.o \
X    obliq.o \
X    parallax.o \
X    pelement.o \
X    plans.o \
X    plot.o \
X    popup.o \
X    precess.o \
X    reduce.o \
X    refract.o \
X    riset.o \
X    riset_c.o \
X    sel_fld.o \
X    sex_dec.o \
X    srch.o \
X    sun.o \
X    time.o \
X    utc_gst.o \
X    version.o \
X    watch.o
X
Xephem:    $(EPHEM)
X    cc -o $@ $(EPHEM) -ltermcap -lm
EOFxEOF
len=`wc -c < Makefile`
if expr $len != 484 > /dev/null
then echo Length of Makefile is $len but it should be 484.
fi
echo x ephem.cfg
sed -e 's/^X//' << 'EOFxEOF' > ephem.cfg
XUT=1:45:0
XUD=4/5/1990
XTZNAME=CST
XTZONE=6:0;0
XLONG=93:42:8
XLAT=44:50:37
XHEIGHT=800
XTEMP=40
XPRES=29.5
XSTPSZ=RTC
XPROPTS=TSMevmjsunpx
XEPOCH=Eod
XNSTEP=1
X
X* comet halley. elements from Duffett-Smith book; mag from 12/86 S&T pg 666
XOBJX=e,1986.109,76.0081,170.011,0.9673,162.2384,58.1540,17.9435,3.66,7.05
X
X* orion, roughly
XOBJX=fixed,6:0:0,0:0:0,1950
X
X* comet austin, 1989c1, as per IAU Circular 4941 and magnitude rumors.
XOBJX=parabolic,4/9.847/1990,58.911,61.504,.34963,75.409,1950,3.8,13.7
EOFxEOF
len=`wc -c < ephem.cfg`
if expr $len != 487 > /dev/null
then echo Length of ephem.cfg is $len but it should be 487.
fi
echo x astro.h
sed -e 's/^X//' << 'EOFxEOF' > astro.h
X#ifndef PI
X#define    PI        3.141592653589793
X#endif
X
X/* conversions among hours (of ra), degrees and radians. */
X#define    degrad(x)    ((x)*PI/180.)
X#define    raddeg(x)    ((x)*180./PI)
X#define    hrdeg(x)    ((x)*15.)
X#define    deghr(x)    ((x)/15.)
X#define    hrrad(x)    degrad(hrdeg(x))
X#define    radhr(x)    deghr(raddeg(x))
X
X/* ratio of from synodic (solar) to sidereal (stellar) rate */
X#define    SIDRATE        .9972695677
X
X/* manifest names for planets.
X * N.B. must cooincide with usage in pelement.c and plans.c.
X */
X#define    MERCURY    0
X#define    VENUS    1
X#define    MARS    2
X#define    JUPITER    3
X#define    SATURN    4
X#define    URANUS    5
X#define    NEPTUNE    6
X#define    PLUTO    7
EOFxEOF
len=`wc -c < astro.h`
if expr $len != 620 > /dev/null
then echo Length of astro.h is $len but it should be 620.
fi
echo x circum.h
sed -e 's/^X//' << 'EOFxEOF' > circum.h
X#define    SPD    (24.0*3600.0)    /* seconds per day */
X
X#define    EOD    (-9786)        /* special epoch flag: use epoch of date */
X#define    RTC    (-1324)        /* special tminc flag: use rt clock */
X#define    NOMAG    (-2134)        /* special s_mag flag: means undefined */
X#define    NOHELIO    (-2314)        /* special s_hlong flag: means it and s_hlat are
X                 * undefined
X                 */
X
X#define    STDHZN        0    /* rise/set times based on nominal conditions */
X#define    ADPHZN        1    /* rise/set times based on exact current " */
X#define    TWILIGHT    2    /* rise/set times for sun 18 degs below hor */
X
X/* info about our local observing circumstances */
Xtypedef struct {
X    double n_mjd;    /* modified Julian date, ie, days since
X             * Jan 0.5 1900 (== 12 noon, Dec 30, 1899), utc.
X             * enough precision to get well better than 1 second.
X             */
X    double n_lat;    /* latitude, >0 north, rads */
X    double n_lng;    /* longitude, >0 east, rads */
X    double n_tz;    /* time zone, hrs behind UTC */
X    double n_temp;    /* atmospheric temp, degrees C */
X    double n_pressure; /* atmospheric pressure, mBar */
X    double n_height;    /* height above sea level, earth radii */
X    double n_epoch;    /* desired precession display epoch as an mjd, or EOD */
X    char n_tznm[4];    /* time zone name; 3 chars or less, always 0 at end */
X} Now;
Xextern double    mjd_day(), mjd_hr();
X
X/* info about where and how we see something in the sky */
Xtypedef struct {
X    double s_ra;    /* ra, rads (precessed to n_epoch) */
X    double s_dec;    /* dec, rads (precessed to n_epoch) */
X    double s_az;    /* azimuth, >0 e of n, rads */
X    double s_alt;    /* altitude above topocentric horizon, rads */
X    double s_sdist;    /* dist from object to sun, au */
X    double s_edist;    /* dist from object to earth, au */
X    double s_elong;    /* angular sep between object and sun, >0 if east */
X    double s_hlong;    /* heliocentric longitude, rads */
X    double s_hlat;    /* heliocentric latitude, rads */
X    double s_size;    /* angular size, arc secs */
X    double s_phase;    /* phase, % */
X    double s_mag;    /* visual magnitude */
X} Sky;
X
X/* flags for riset_cir() status */
X#define    RS_NORISE    0x001    /* object does not rise as such today */
X#define    RS_2RISES    0x002    /* object rises more than once today */
X#define    RS_NOSET    0x004    /* object does not set as such today */
X#define    RS_2SETS    0x008    /* object sets more than once today */
X#define    RS_CIRCUMPOLAR    0x010    /* object stays up all day today */
X#define    RS_2TRANS    0x020    /* transits twice in one day */
X#define    RS_NEVERUP    0x040    /* object never rises today */
X#define    RS_NOTRANS    0x080    /* doesn't transit today */
X#define    RS_ERROR    0x100    /* can't figure out times... */
X
X/* shorthands for fields a Now pointer, np */
X#define mjd    np->n_mjd
X#define lat    np->n_lat
X#define lng    np->n_lng
X#define tz    np->n_tz
X#define temp    np->n_temp
X#define pressure np->n_pressure
X#define height    np->n_height
X#define epoch    np->n_epoch
X#define tznm    np->n_tznm
EOFxEOF
len=`wc -c < circum.h`
if expr $len != 2795 > /dev/null
then echo Length of circum.h is $len but it should be 2795.
fi
echo x screen.h
sed -e 's/^X//' << 'EOFxEOF' > screen.h
X/* screen layout details
X *
X * it looks better if the fields are drawn in some nice order so it you
X * rearrange the fields, check the menu printing functions.
X */
X
X/* size of screen */
X#define    NR    24
X#define    NC    80
X
X#define    ASPECT    (4./3.)    /* screen width to height dimensions ratio */
X
X#define    GAP    6    /* gap between field name and value */
X
X#define    COL1        1
X#define    COL2        27
X#define    COL3        44
X#define    COL4        61    /* calendar */
X
X#define    R_PROMPT    1    /* prompt row */
X#define    C_PROMPT    COL1
X
X#define    R_NEWCIR    2
X#define    C_NEWCIR    ((NC-17)/2) /* 17 is length of the message */
X
X#define    R_TOP        3    /* first row of top menu items */
X
X#define    R_TZN    (R_TOP+0)
X#define    C_TZN    COL1
X#define    R_LT    R_TZN
X#define    C_LT    (C_TZN+GAP-2)
X#define    R_LD    R_TZN
X#define    C_LD    (C_TZN+13)
X
X#define    R_UT    (R_TOP+1)
X#define    C_UT    COL1
X#define    C_UTV    (C_UT+GAP-2)
X#define    R_UD    R_UT
X#define    C_UD    (C_UT+13)
X
X#define    R_JD    (R_TOP+2)
X#define    C_JD    COL1
X#define    C_JDV    (C_JD+GAP+3)
X
X#define    R_LST    (R_TOP)
X#define    C_LST    COL2
X#define    C_LSTV    (C_LST+GAP)
X
X#define    R_LAT    (R_TOP+0)
X#define    C_LAT    COL3
X#define    C_LATV    (C_LAT+4)
X
X#define    R_DAWN    (R_TOP+2)
X#define    C_DAWN    COL2
X#define    C_DAWNV    (C_DAWN+GAP+3)
X
X#define    R_STPSZ    (R_TOP+7)
X#define    C_STPSZ    COL2
X#define    C_STPSZV (C_STPSZ+GAP-1)
X
X#define    R_HEIGHT (R_TOP+2)
X#define    C_HEIGHT COL3
X#define    C_HEIGHTV (C_HEIGHT+GAP)
X
X#define    R_PRES    (R_TOP+4)
X#define    C_PRES    COL3
X#define    C_PRESV    (C_PRES+GAP)
X
X#define    R_WATCH    (R_TOP+4)
X#define    C_WATCH    COL1
X
X#define    R_SRCH    (R_TOP+5)
X#define    C_SRCH    COL1
X#define    C_SRCHV    (C_SRCH+16)
X
X#define    R_PLOT    (R_TOP+6)
X#define    C_PLOT    (COL1)
X#define    C_PLOTV (C_PLOT+20)
X
X#define    R_ALTM    (R_TOP+7)
X#define    C_ALTM    COL1
X#define    C_ALTMV    (C_ALTM+10)
X
X#define    R_TZONE    (R_TOP+5)
X#define    C_TZONE    COL3
X#define    C_TZONEV (C_TZONE+GAP-1)
X
X#define    R_LONG    (R_TOP+1)
X#define    C_LONG    COL3
X#define    C_LONGV    (C_LONG+4)
X
X#define    R_DUSK    (R_TOP+3)
X#define    C_DUSK    COL2
X#define    C_DUSKV    (C_DUSK+GAP+3)
X
X#define    R_NSTEP (R_TOP+6)
X#define    C_NSTEP    COL2
X#define    C_NSTEPV (C_NSTEP+GAP)
X
X#define    R_TEMP    (R_TOP+3)
X#define    C_TEMP    COL3
X#define    C_TEMPV    (C_TEMP+GAP)
X
X#define    R_EPOCH        (R_TOP+6)
X#define    C_EPOCH        COL3
X#define    C_EPOCHV    (C_EPOCH+GAP)
X
X#define    R_MNUDEP    (R_TOP+6)
X#define    C_MNUDEP    COL3
X#define    C_MNUDEPV    (C_EPOCH+GAP)
X
X#define    R_LON    (R_TOP+4)
X#define    C_LON    COL2
X#define    C_LONV    (C_LON+GAP+3)
X
X#define    R_CAL    R_TOP
X#define    C_CAL   COL4
X
X/* menu 1 info table */
X#define    R_PLANTAB    (R_TOP+9)
X#define    R_SUN        (R_PLANTAB+2)
X#define    R_MOON        (R_PLANTAB+3)
X#define    R_MERCURY    (R_PLANTAB+4)
X#define    R_VENUS        (R_PLANTAB+5)
X#define    R_MARS        (R_PLANTAB+6)
X#define    R_JUPITER    (R_PLANTAB+7)
X#define    R_SATURN    (R_PLANTAB+8)
X#define    R_URANUS    (R_PLANTAB+9)
X#define    R_NEPTUNE    (R_PLANTAB+10)
X#define    R_PLUTO        (R_PLANTAB+11)
X#define    R_OBJX        (R_PLANTAB+12)
X#define    C_OBJ        1
X#define    C_RA        4
X#define    C_DEC        12
X#define    C_AZ        19
X#define    C_ALT        26
X#define    C_HLONG        33
X#define    C_HLAT        40
X#define    C_EDIST        47
X#define C_SDIST     54
X#define    C_ELONG        61
X#define    C_SIZE        68
X#define    C_MAG        73
X#define    C_PHASE        78
X
X/* menu 2 screen items */
X#define    C_RISETM    10
X#define    C_RISEAZ    18
X#define    C_TRANSTM    31
X#define    C_TRANSALT    39
X#define    C_SETTM        52
X#define    C_SETAZ        60
X#define    C_TUP        72
X
X/* menu 3 items */
X#define    C_SUN        4
X#define    C_MOON        11
X#define    C_MERCURY    18
X#define    C_VENUS        25
X#define    C_MARS        32
X#define    C_JUPITER    39
X#define    C_SATURN    46
X#define    C_URANUS    53
X#define    C_NEPTUNE    60
X#define    C_PLUTO        67
X#define    C_OBJX        74
X
X#define    PW    (NC-C_PROMPT+1)    /* total prompt line width */
X
X/* macros to pack a row/col and menu selection flags all into an int.
X * (use this rather than a structure because we can compare them so easily.
X * could use bit fields and a union, but then can't init them or use switch.)
X * bit field defs: [15..14]=menu  [13..12]=flags  [11..7]=row  [6..0]=column.
X * see sel_fld.c.
X * F_MNUX also used in main to manage which bottom menu is up.
X */
X#define    F_CHG        (1<<12)    /* field may be picked for changing */
X#define    F_PLT        (1<<13)    /* field may be picked for plotting */
X#define    F_MMNU        (0<<14)    /* field is on main menu */
X#define    F_MNU1        (1<<14)    /* field is on menu 1 */
X#define    F_MNU2         (2<<14)    /* field is on menu 2 */
X#define    F_MNU3        (3<<14)    /* field is on menu 3 */
X#define    rcfpack(r,c,f)    ((f) | ((r) << 7) | (c))
X#define    unpackr(p)    (((p) >> 7) & 0x1f)
X#define    unpackc(p)    ((p) & 0x7f)
X#define    unpackrc(p)    ((p) & 0xfff)
X#define    tstpackf(p,f)    (((p) & ((f)&0x3000)) && \
X            (((p)&0xc000) == ((f)&0xc000) || ((p)&0xc000) == 0))
X
X/* additions to the planet defines from astro.h.
X * must not conflict, and must fit in range 0..15.
X */
X#define    SUN    (PLUTO+1)
X#define    MOON    (PLUTO+2)
X#define    OBJX    (PLUTO+3)    /* the user-defined object */
X#define    NOBJ    (OBJX+1)    /* total number of objects */
X
X#define    cntrl(x)    ((x) & 037)
X#define    QUIT        cntrl('d')    /* char to exit program */
X#define    HELP        '?'        /* char to give help message */
X#define    REDRAW        cntrl('l')    /* char to redraw (like vi) */
X#define    VERSION        cntrl('v')    /* char to display version number */
X#define    END        'q'        /* char to quit current mode */
EOFxEOF
len=`wc -c < screen.h`
if expr $len != 4950 > /dev/null
then echo Length of screen.h is $len but it should be 4950.
fi
echo x aa_hadec.c
sed -e 's/^X//' << 'EOFxEOF' > aa_hadec.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
X * azimuth (angle round to the east from north+, radians),
X * return hour angle (radians), ha, and declination (radians), dec.
X */
Xaa_hadec (lat, alt, az, ha, dec)
Xdouble lat;
Xdouble alt, az;
Xdouble *ha, *dec;
X{
X    aaha_aux (lat, az, alt, ha, dec);
X}
X
X/* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
X * (radians), dec,
X * return altitude (up+, radians), alt, and
X * azimuth (angle round to the east from north+, radians),
X */
Xhadec_aa (lat, ha, dec, alt, az)
Xdouble lat;
Xdouble ha, dec;
Xdouble *alt, *az;
X{
X    aaha_aux (lat, ha, dec, az, alt);
X}
X
X/* the actual formula is the same for both transformation directions so
X * do it here once for each way.
X * N.B. all arguments are in radians.
X */
Xstatic
Xaaha_aux (lat, x, y, p, q)
Xdouble lat;
Xdouble x, y;
Xdouble *p, *q;
X{
X    static double lastlat = -1000.;
X    static double sinlastlat, coslastlat;
X    double sy, cy;
X    double sx, cx;
X    double sq, cq;
X    double a;
X    double cp;
X
X    /* latitude doesn't change much, so try to reuse the sin and cos evals.
X     */
X    if (lat != lastlat) {
X        sinlastlat = sin (lat);
X        coslastlat = cos (lat);
X        lastlat = lat;
X    }
X
X    sy = sin (y);
X    cy = cos (y);
X    sx = sin (x);
X    cx = cos (x);
X
X/* define GOODATAN2 if atan2 returns full range -PI through +PI.
X */
X#ifdef GOODATAN2
X    *q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
X    *p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
X#else
X#define    EPS    (1e-20)
X    sq = (sy*sinlastlat) + (cy*coslastlat*cx);
X    *q = asin (sq);
X    cq = cos (*q);
X    a = coslastlat*cq;
X    if (a > -EPS && a < EPS)
X        a = a < 0 ? -EPS : EPS; /* avoid / 0 */
X    cp = (sy - (sinlastlat*sq))/a;
X    if (cp >= 1.0)    /* the /a can be slightly > 1 */
X        *p = 0.0;
X    else if (cp <= -1.0)
X        *p = PI;
X    else
X        *p = acos ((sy - (sinlastlat*sq))/a);
X    if (sx>0) *p = 2.0*PI - *p;
X#endif
X}
EOFxEOF
len=`wc -c < aa_hadec.c`
if expr $len != 1922 > /dev/null
then echo Length of aa_hadec.c is $len but it should be 1922.
fi
echo x altmenus.c
sed -e 's/^X//' << 'EOFxEOF' > altmenus.c
X/* printing routines for the three alternative bottom half menus,
X * "menu1", "menu2" and "menu3".
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h"
X
Xstatic int altmenu = F_MNU1;    /* which alternate menu is up; one of F_MNUi */
Xstatic int alt2_stdhzn;    /* whether to use STDHZN (aot ADPHZN) horizon algthm  */
Xstatic int alt3_geoc;    /* whether to use geocentric (aot topocentric) vantage*/
X
X/* table of screen rows given a body #define from astro/h or screen.h */
Xstatic short bodyrow[NOBJ] = {
X    R_MERCURY, R_VENUS, R_MARS, R_JUPITER, R_SATURN,
X    R_URANUS, R_NEPTUNE, R_PLUTO, R_SUN, R_MOON, R_OBJX
X};
X/* table of screen cols for third menu format, given body #define ... */
Xstatic short bodycol[NOBJ] = {
X    C_MERCURY, C_VENUS, C_MARS, C_JUPITER, C_SATURN,
X    C_URANUS, C_NEPTUNE, C_PLUTO, C_SUN, C_MOON, C_OBJX
X};
X
X/* let op decide which alternate menu should be up,
X * including any menu-specific setup they might require.
X * return 0 if things changed to require updating the alt menu; else -1.
X */
Xaltmenu_setup()
X{
X    static char *flds[5] = {
X        "Data menu", "Rise/Set menu", "Separations menu"
X    };
X    int newmenu = altmenu, newhzn = alt2_stdhzn, newgeoc = alt3_geoc;
X    int new;
X    int fn = altmenu == F_MNU1 ? 0 : altmenu == F_MNU2 ? 1 : 2;
X
X    ask:
X    flds[3]= newhzn ? "Standard hzn" : "Adaptive hzn";
X    flds[4]= newgeoc? "Geocentric" : "Topocentric";
X
X    switch (popup (flds, fn, 5)) {
X    case 0: newmenu = F_MNU1; break;
X    case 1: newmenu = F_MNU2; break;
X    case 2: newmenu = F_MNU3; break;
X    case 3: newhzn ^= 1; fn = 3; goto ask;
X    case 4: newgeoc ^= 1; fn = 4; goto ask;
X    default: return (-1);
X    }
X
X    new = 0;
X    if (newmenu != altmenu) {
X        altmenu = newmenu;
X        new++;
X    }
X    if (newhzn != alt2_stdhzn) {
X        alt2_stdhzn = newhzn;
X        if (newmenu == F_MNU2)
X        new++;
X    }
X    if (newgeoc != alt3_geoc) {
X        alt3_geoc = newgeoc;
X        if (newmenu == F_MNU3)
X        new++;
X    }
X    return (new ? 0 : -1);
X}
X
X/* erase the info for the given planet */
Xalt_nobody (p)
Xint p;
X{
X    f_eol (bodyrow[p], C_RA);
X}
X
Xalt_body (b, force, np)
Xint b;        /* which body, ala astro.h and screen.h defines */
Xint force;    /* if !0 then draw for sure, else just if changed since last */
XNow *np;
X{
X    switch (altmenu) {
X    case F_MNU1: alt1_body (b, force, np); break;
X    case F_MNU2: alt2_body (b, force, np); break;
X    case F_MNU3: alt3_body (b, force, np); break;
X    }
X}
X
X/* draw the labels for the current alternate menu format */
Xalt_labels ()
X{
X    switch (altmenu) {
X    case F_MNU1: alt1_labels (); break;
X    case F_MNU2: alt2_labels (); break;
X    case F_MNU3: alt3_labels (); break;
X    }
X}
X
X/* erase the labels for the current alternate menu format */
Xalt_nolabels ()
X{
X    int i;
X
X    f_string (R_ALTM, C_ALTMV, "             ");
X    for (i = R_PLANTAB; i < R_SUN; i++)
X        f_eol (i, C_RA);
X}
X
Xalt_menumask()
X{
X    return (altmenu);
X}
X
X/* handy function to return the next planet in the order in which they are
X * displayed in the lower half of the screen.
X * input is a given planet, return is the next planet.
X * if input is not legal, then first planet is returned; when input is the
X * last planet, then -1 is returned.
X * typical usage is something like:
X *   for (p = nxtbody(-1); p != -1; p = nxtbody(p))
X */
Xnxtbody(c)
Xint c;
X{
X    static short nxtpl[NOBJ] = {
X        VENUS, MARS, JUPITER, SATURN, URANUS,
X        NEPTUNE, PLUTO, OBJX, MOON, MERCURY, -1
X    };
X
X    if (c < MERCURY || c >= NOBJ)
X        return (SUN);
X    else
X        return (nxtpl[c]);
X}
X
Xstatic
Xalt1_labels()
X{
X    f_string (R_ALTM, C_ALTMV, "  Planet Data");
X
X    f_string (R_PLANTAB,    C_RA,    "R.A.");
X    f_string (R_PLANTAB+1,    C_RA,    "Hr:Mn.d");
X    f_string (R_PLANTAB,    C_DEC,    "Dec");
X    f_string (R_PLANTAB+1,    C_DEC,    "Deg:Mn");
X    f_string (R_PLANTAB,    C_AZ,    "Az");
X    f_string (R_PLANTAB+1,    C_AZ,    "Deg E");
X    f_string (R_PLANTAB,    C_ALT,    "Alt");
X    f_string (R_PLANTAB+1,    C_ALT,    "Deg Up");
X    f_string (R_PLANTAB,    C_HLONG,"Helio");
X    f_string (R_PLANTAB+1,    C_HLONG,"Long");
X    f_string (R_PLANTAB,    C_HLAT,    "Helio");
X    f_string (R_PLANTAB+1,    C_HLAT,    "Lat");
X    f_string (R_PLANTAB,    C_EDIST,"Ea Dst");
X    f_string (R_PLANTAB+1,    C_EDIST,"AU(mi)");
X    f_string (R_PLANTAB,    C_SDIST,"Sn Dst");
X    f_string (R_PLANTAB+1,    C_SDIST,"AU");
X    f_string (R_PLANTAB,    C_ELONG,"Elong");
X    f_string (R_PLANTAB+1,    C_ELONG,"Deg E");
X    f_string (R_PLANTAB,    C_SIZE,    "Size");
X    f_string (R_PLANTAB+1,    C_SIZE,    "ArcS");
X    f_string (R_PLANTAB,    C_MAG,    "VMag");
X    f_string (R_PLANTAB,    C_PHASE,"Phs");
X    f_char (R_PLANTAB+1,    C_PHASE,'%');
X}
X
Xstatic
Xalt2_labels()
X{
X    f_string (R_ALTM, C_ALTMV, "Rise/Set Info");
X
X    f_string (R_PLANTAB,    C_RISETM+5,    "Rise");
X    f_string (R_PLANTAB+1,    C_RISETM+1,    "Time");
X    f_string (R_PLANTAB+1,    C_RISEAZ+2,    "Az");
X    f_string (R_PLANTAB,    C_TRANSTM+3,    "Transit");
X    f_string (R_PLANTAB+1,    C_TRANSTM+1,    "Time");
X    f_string (R_PLANTAB+1,    C_TRANSALT+2,    "Alt");
X    f_string (R_PLANTAB,    C_SETTM+5,    "Set");
X    f_string (R_PLANTAB+1,    C_SETTM+1,    "Time");
X    f_string (R_PLANTAB+1,    C_SETAZ+2,    "Az");
X    f_string (R_PLANTAB,    C_TUP,        "Hrs Up");
X}
X
Xstatic
Xalt3_labels()
X{
X    f_string (R_ALTM, C_ALTMV, "  Separations");
X
X    f_string (R_PLANTAB,    C_SUN+2,    "Sun");
X    f_string (R_PLANTAB,    C_MOON+1,    "Moon");
X    f_string (R_PLANTAB,    C_MERCURY+1,    "Merc");
X    f_string (R_PLANTAB,    C_VENUS+1,    "Venus");
X    f_string (R_PLANTAB,    C_MARS+1,    "Mars");
X    f_string (R_PLANTAB,    C_JUPITER+2,    "Jup");
X    f_string (R_PLANTAB,    C_SATURN,    "Saturn");
X    f_string (R_PLANTAB,    C_URANUS,    "Uranus");
X    f_string (R_PLANTAB,    C_NEPTUNE+2,    "Nep");
X    f_string (R_PLANTAB,    C_PLUTO+1,    "Pluto");
X    f_string (R_PLANTAB,    C_OBJX+3,    "X");
X}
X
X/* print body info in first menu format */
Xstatic
Xalt1_body (p, force, np)
Xint p;        /* which body, as in astro.h/screen.h defines */
Xint force;    /* whether to print for sure or only if things have changed */
XNow *np;
X{
X    Sky sky;
X    double as = plot_ison() || srch_ison() ? 0.0 : 60.0;
X    int row = bodyrow[p];
X
X    if (body_cir (p, as, np, &sky) || force) {
X        f_ra (row, C_RA, sky.s_ra);
X        f_angle (row, C_DEC, sky.s_dec);
X        if (p != MOON && sky.s_hlong != NOHELIO) {
X        f_angle (row, C_HLONG, sky.s_hlong);
X        if (p != SUN)
X            f_angle (row, C_HLAT, sky.s_hlat);
X        }
X
X        if (p == MOON) {
X        /* distance is on km, show in miles */
X        f_double (R_MOON, C_EDIST, "%6.0f", sky.s_edist/1.609344);
X        } else if (sky.s_edist > 0.0) {
X        /* show distance in au */
X        f_double (row, C_EDIST,(sky.s_edist>=10.0)?"%6.3f":"%6.4f",
X                                sky.s_edist);
X        }
X        if (sky.s_sdist > 0.0)
X        f_double (row, C_SDIST, (sky.s_sdist>=10.0)?"%6.3f":"%6.4f",
X                                sky.s_sdist);
X        if (p != SUN)
X        f_double (row, C_ELONG, "%6.1f", sky.s_elong);
X        if (p != OBJX)
X        f_double (row, C_SIZE, (p==MOON||p==SUN)?"%4.0f":"%4.1f",
X                                sky.s_size);
X        if (sky.s_mag != NOMAG)
X        f_double (row, C_MAG, (p==MOON||p==SUN)?"%4.0f":"%4.1f",
X                                    sky.s_mag);
X        if (sky.s_sdist > 0.0) {
X        char buf[10];
X        /* some terminals scroll when write a char in low-right corner.
X         * TODO: is there a nicer way to handle this maybe?
X         */
X        int col = row == NR ? C_PHASE - 1 : C_PHASE;
X        /* would just do this if Turbo-C 2.0 "%?.0f" worked:
X         * f_double (row, col, "%3.0f", sky.s_phase);
X         */
X        sprintf (buf, "%3d", (int)(sky.s_phase+0.5));
X        f_string (row, col, buf);
X        }
X    }
X
X    f_angle (row, C_AZ, sky.s_az);
X    f_angle (row, C_ALT, sky.s_alt);
X}
X
X/* print body info in the second menu format */
Xstatic
Xalt2_body (p, force, np)
Xint p;        /* which body, as in astro.h/screen.h defines */
Xint force;    /* whether to print for sure or only if things have changed */
XNow *np;
X{
X    double ltr, lts, ltt, azr, azs, altt;
X    int row = bodyrow[p];
X    int status;
X    double tmp;
X    int today_tup = 0;
X
X    if (!riset_cir (p, np, alt2_stdhzn?STDHZN:ADPHZN, <r, <s, <t,
X                    &azr, &azs, &altt, &status) && !force)
X        return;
X
X    alt_nobody (p);
X
X    if (status & RS_ERROR) {
X        /* can not find where body is! */
X        f_string (row, C_RISETM, "?Error?");
X        return;
X    }
X    if (status & RS_CIRCUMPOLAR) {
X        /* body is up all day */
X        f_string (row, C_RISETM, "Circumpolar");
X        if (status & RS_NOTRANS)
X        f_string (row, C_TRANSTM, "No transit");
X        else {
X        f_mtime (row, C_TRANSTM, ltt);
X        if (status & RS_2TRANS)
X            f_char (row, C_TRANSTM+5, '+');
X        f_angle (row, C_TRANSALT, altt);
X        }
X        f_string (row, C_TUP, "24:00"); /*f_mtime() changes to 0:00 */
X        return;
X    }
X    if (status & RS_NEVERUP) {
X        /* body never up at all today */
X        f_string (row, C_RISETM, "Never up");
X        f_mtime (row, C_TUP, 0.0);
X        return;
X    }
X
X    if (status & RS_NORISE) {
X        /* object does not rise as such today */
X        f_string (row, C_RISETM, "Never rises");
X        ltr = 0.0; /* for TUP */
X        today_tup = 1;
X    } else {
X        f_mtime (row, C_RISETM, ltr);
X        if (status & RS_2RISES) {
X        /* object rises more than once today */
X        f_char (row, C_RISETM+5, '+');
X        }
X        f_angle (row, C_RISEAZ, azr);
X    }
X
X    if (status & RS_NOTRANS)
X        f_string (row, C_TRANSTM, "No transit");
X    else {
X        f_mtime (row, C_TRANSTM, ltt);
X        if (status & RS_2TRANS)
X        f_char (row, C_TRANSTM+5, '+');
X        f_angle (row, C_TRANSALT, altt);
X    }
X
X    if (status & RS_NOSET) {
X        /* object does not set as such today */
X        f_string (row, C_SETTM, "Never sets");
X        lts = 24.0;    /* for TUP */
X        today_tup = 1;
X    } else {
X        f_mtime (row, C_SETTM, lts);
X        if (status & RS_2SETS)
X        f_char (row, C_SETTM+5, '+');
X        f_angle (row, C_SETAZ, azs);
X    }
X
X    tmp = lts - ltr;
X    if (tmp < 0)
X        tmp = 24.0 + tmp;
X    f_mtime (row, C_TUP, tmp);
X    if (today_tup)
X        f_char (row, C_TUP+5, '+');
X}
X
X/* print body info in third menu format. this may be either the geocentric
X *   or topocentric angular separation between object p and each of the others.
X *   the latter, of course, includes effects of refraction and so can change
X *   quite rapidly near the time of each planets rise or set.
X * for now, we don't save old values so we always redo everything and ignore
X *  the "force" argument. this isn't that bad since body_cir() has memory and
X *   will avoid most computations as we hit them again in the lower triangle.
X */
X/*ARGSUSED*/
Xstatic
Xalt3_body (p, force, np)
Xint p;        /* which body, as in astro.h/screen.h defines */
Xint force;    /* whether to print for sure or only if things have changed */
XNow *np;
X{
X    int row = bodyrow[p];
X    Sky skyp, skyq;
X    double spy, cpy, px, *qx, *qy;
X    int wantx = objx_ison();
X    double as = plot_ison() || srch_ison() ? 0.0 : 60.0;
X    int q;
X
X    (void) body_cir (p, as, np, &skyp);
X    if (alt3_geoc) {
X        /* use ra for "x", dec for "y". */
X        spy = sin (skyp.s_dec);
X        cpy = cos (skyp.s_dec);
X        px = skyp.s_ra;
X        qx = &skyq.s_ra;
X        qy = &skyq.s_dec;
X    } else {
X        /* use azimuth for "x", altitude for "y". */
X        spy = sin (skyp.s_alt);
X        cpy = cos (skyp.s_alt);
X        px = skyp.s_az;
X        qx = &skyq.s_az;
X        qy = &skyq.s_alt;
X    }
X    for (q = nxtbody(-1); q != -1; q = nxtbody(q))
X        if (q != p && (q != OBJX || wantx)) {
X        double csep;
X        (void) body_cir (q, as, np, &skyq);
X        csep = spy*sin(*qy) + cpy*cos(*qy)*cos(px-*qx);
X        f_angle (row, bodycol[q], acos(csep));
X        }
X}
EOFxEOF
len=`wc -c < altmenus.c`
if expr $len != 10912 > /dev/null
then echo Length of altmenus.c is $len but it should be 10912.
fi
echo x anomaly.c
sed -e 's/^X//' << 'EOFxEOF' > anomaly.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define    TWOPI    (2*PI)
X
X/* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
X * find the true anomaly, *nu, and the eccentric anomaly, *ea.
X * all angles in radians.
X */
Xanomaly (ma, s, nu, ea)
Xdouble ma, s;
Xdouble *nu, *ea;
X{
X    double m, dla, fea;
X
X    m = ma-TWOPI*(long)(ma/TWOPI);
X    fea = m;
X    while (1) {
X        dla = fea-(s*sin(fea))-m;
X        if (fabs(dla)<1e-6)
X        break;
X        dla /= 1-(s*cos(fea));
X        fea -= dla;
X    }
X
X    *nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
X    *ea = fea;
X}
EOFxEOF
len=`wc -c < anomaly.c`
if expr $len != 557 > /dev/null
then echo Length of anomaly.c is $len but it should be 557.
fi
echo x cal_mjd.c
sed -e 's/^X//' << 'EOFxEOF' > cal_mjd.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given a date in months, mn, days, dy, years, yr,
X * return the modified Julian date (number of days elapsed since 1900 jan 0.5),
X * *mjd.
X */
Xcal_mjd (mn, dy, yr, mjd)
Xint mn, yr;
Xdouble dy;
Xdouble *mjd;
X{
X    int b, d, m, y;
X    long c;
X
X    m = mn;
X    y = (yr < 0) ? yr + 1 : yr;
X    if (mn < 3) {
X        m += 12;
X        y -= 1;
X    }
X
X    if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15))
X        b = 0;
X    else {
X        int a;
X        a = y/100;
X        b = 2 - a + a/4;
X    }
X
X    if (y < 0)
X        c = (long)((365.25*y) - 0.75) - 694025L;
X    else
X        c = (long)(365.25*y) - 694025L;
X
X    d = 30.6001*(m+1);
X
X    *mjd = b + c + d + dy - 0.5;
X}
X
X/* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
X * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
X */
Xmjd_cal (mjd, mn, dy, yr)
Xdouble mjd;
Xint *mn, *yr;
Xdouble *dy;
X{
X    double d, f;
X    double i, a, b, ce, g;
X
X    d = mjd + 0.5;
X    i = floor(d);
X    f = d-i;
X    if (f == 1) {
X        f = 0;
X        i += 1;
X    }
X
X    if (i > -115860.0) {
X        a = floor((i/36524.25)+.9983573)+14;
X        i += 1 + a - floor(a/4.0);
X    }
X
X    b = floor((i/365.25)+.802601);
X    ce = i - floor((365.25*b)+.750001)+416;
X    g = floor(ce/30.6001);
X    *mn = g - 1;
X    *dy = ce - floor(30.6001*g)+f;
X    *yr = b + 1899;
X
X    if (g > 13.5)
X        *mn = g - 13;
X    if (*mn < 2.5)
X        *yr = b + 1900;
X    if (*yr < 1)
X        *yr -= 1;
X}
X
X/* given an mjd, set *dow to 0..6 according to which dayof the week it falls
X * on (0=sunday) or set it to -1 if can't figure it out.
X */
Xmjd_dow (mjd, dow)
Xdouble mjd;
Xint *dow;
X{
X    /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
X     * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
X     * year algorithm). however, Great Britian and the colonies did not
X     * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
X     * due to additional accumulated error). leap years before 1752 thus
X     * can not easily be accounted for from the cal_mjd() number...
X     */
X    if (mjd < -53798.5) {
X        /* pre sept 14, 1752 too hard to correct */
X        *dow = -1;
X        return;
X    }
X    *dow = ((long)floor(mjd-.5) + 1) % 7;/* 1/1/1900 (mjd 0.5) is a Monday*/
X    if (*dow < 0)
X        *dow += 7;
X}
X
X/* given a mjd, return the the number of days in the month.  */
Xmjd_dpm (mjd, ndays)
Xdouble mjd;
Xint *ndays;
X{
X    static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
X    int m, y;
X    double d;
X
X    mjd_cal (mjd, &m, &d, &y);
X    *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
X}
X
X
X/* given a mjd, return the year as a double. */
Xmjd_year (mjd, yr)
Xdouble mjd;
Xdouble *yr;
X{
X    int m, y;
X    double d;
X    double e0, e1;    /* mjd of start of this year, start of next year */
X
X    mjd_cal (mjd, &m, &d, &y);
X    cal_mjd (1, 1.0, y, &e0);
X    cal_mjd (1, 1.0, y+1, &e1);
X    *yr = y + (mjd - e0)/(e1 - e0);
X}
X
X/* given a decimal year, return mjd */
Xyear_mjd (y, mjd)
Xdouble y;
Xdouble *mjd;
X{
X    double e0, e1;    /* mjd of start of this year, start of next year */
X    int yf = floor (y);
X
X    cal_mjd (1, 1.0, yf, &e0);
X    cal_mjd (1, 1.0, yf+1, &e1);
X    *mjd = e0 + (y - yf)*(e1-e0);
X}
EOFxEOF
len=`wc -c < cal_mjd.c`
if expr $len != 3067 > /dev/null
then echo Length of cal_mjd.c is $len but it should be 3067.
fi
echo x circum.c
sed -e 's/^X//' << 'EOFxEOF' > circum.c
X/* fill in a Sky struct with all we know about each object.
X *(object-x is in objx.c)
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h"    /* just for SUN and MOON */
X
X/* find body p's circumstances now.
X * to save some time the caller may specify a desired accuracy, in arc seconds.
X * if, based on its mean motion, it would not have moved this much since the
X * last time we were called we only recompute altitude and azimuth and avoid
X * recomputing the planet's heliocentric position. use 0.0 for best possible.
X * we always recompute objx's position regardless.
X * return 0 if only alt/az changes, else 1 if all other stuff updated too.
X * TODO: beware of fast retrograde motions.
X */
Xbody_cir (p, as, np, sp)
Xint p;
Xdouble as;
XNow *np;
XSky *sp;
X{
X    typedef struct {
X        double l_dpas;    /* mean days per arc second */
X        Now l_now;        /* when l_sky was found */
X        double l_ra, l_dec;    /* the eod, ie, unprecessed, ra/dec values */
X        Sky l_sky;
X    } Last;
X    /* must be in same order as the astro.h object #define's */
X    static Last last[8] =
X        {{.000068},{.00017},{.00053},{.0034},{.0092},{.027},{.046},{.069}};
X    Last objxlast;
X    double lst, alt, az;
X    double ehp, ha, dec;    /* ehp: angular dia of earth from body */
X    Last *lp;
X    int new;
X
X    switch (p) {
X    case SUN: return (sun_cir (as, np, sp));
X    case MOON: return (moon_cir (as, np, sp));
X    case OBJX: lp = &objxlast; break;
X    default: lp = last + p; break;
X    }
X
X    /* if less than l_every days from last time for this planet
X     * just redo alt/az.
X     * always redo object x.
X     */
X    if (p != OBJX && same_cir (np, &lp->l_now)
X              && about_now (np, &lp->l_now, as*lp->l_dpas)) {
X        *sp = lp->l_sky;
X        new = 0;
X    } else {
X        double lpd0, psi0;    /* heliocentric ecliptic long and lat */
X        double rp0;        /* dist from sun */
X        double rho0;    /* dist from earth */
X        double lam, bet;    /* geocentric ecliptic long and lat */
X        double dia, mag;    /* angular diameter at 1 AU and magnitude */
X        double lsn, rsn;    /* true geoc lng of sun, dist from sn to earth*/
X        double el;    /* elongation */
X        double f;    /* phase from earth */
X
X        lp->l_now = *np;
X        sunpos (mjd, &lsn, &rsn);
X        if (p == OBJX)
X        objx_cir(mjd, &lpd0, &psi0, &rp0, &rho0, &lam, &bet,&sp->s_mag);
X        else {
X        double deps, dpsi;
X        double a;
X        plans(mjd, p, &lpd0, &psi0, &rp0, &rho0, &lam, &bet, &dia,&mag);
X        nutation (mjd, &deps, &dpsi);    /* correct for nutation */
X        lam += dpsi;
X        a = lsn-lam;            /* and 20.4" aberation */
X        lam -= degrad(20.4/3600)*cos(a)/cos(bet);
X        bet -= degrad(20.4/3600)*sin(a)*sin(bet);
X        }
X
X        ecl_eq (mjd, bet, lam, &lp->l_ra, &lp->l_dec);
X
X        sp->s_ra = lp->l_ra;
X        sp->s_dec = lp->l_dec;
X        if (epoch != EOD)
X        precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
X        sp->s_edist = rho0;
X        sp->s_sdist = rp0;
X        elongation (lam, bet, lsn, &el);
X        el = raddeg(el);
X        sp->s_elong = el;
X        sp->s_size = (p != OBJX) ? dia/rho0 : 0.0;
X        f = 0.5*(1+cos(lam-lpd0));
X        sp->s_phase = f*100.0; /* percent */
X        /* objx_cir already returned the apparent magnitude */
X        if (p != OBJX)
X        sp->s_mag = mag + 5.0*log(rp0*rho0/sqrt(f))/log(10.0);
X        sp->s_hlong = lpd0;
X        sp->s_hlat = psi0;
X        new = 1;
X    }
X
X    /* alt, az; correct for parallax and refraction; use eod ra/dec */
X    now_lst (np, &lst);
X    ha = hrrad(lst) - lp->l_ra;
X    if (sp->s_edist > 0.0) {
X        ehp = (2.0*6378.0/146.0e6) / sp->s_edist;
X        ta_par (ha, lp->l_dec, lat, height, ehp, &ha, &dec);
X    } else
X        dec = lp->l_dec;
X    hadec_aa (lat, ha, dec, &alt, &az);
X    refract (pressure, temp, alt, &alt);
X    sp->s_alt = alt;
X    sp->s_az = az;
X    lp->l_sky = *sp;
X    return (new);
X}
X
X/* find local times when sun is 18 degrees below horizon.
X * return 0 if just returned same stuff as previous call, else 1 if new.
X */
Xtwilight_cir (np, dawn, dusk, status)
XNow *np;
Xdouble *dawn, *dusk;
Xint *status;
X{
X    static Now last_now;
X    static double last_dawn, last_dusk;
X    static int last_status;
X    int new;
X
X    if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
X        *dawn = last_dawn;
X        *dusk = last_dusk;
X        *status = last_status;
X        new = 0;
X    } else {
X        double x;
X        (void) riset_cir (SUN,np,TWILIGHT,dawn,dusk,&x,&x,&x,&x,status);
X        last_dawn = *dawn;
X        last_dusk = *dusk;
X        last_status = *status;
X        last_now = *np;
X        new = 1;
X    }
X    return (new);
X}
X
X/* find sun's circumstances now.
X * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
X * return 0 if only alt/az changes, else 1 if all other stuff updated too.
X */
Xsun_cir (as, np, sp)
Xdouble as;
XNow *np;
XSky *sp;
X{
X    static Sky last_sky;
X    static Now last_now;
X    static double last_ra, last_dec;    /* unprecessed ra/dec */
X    double lst, alt, az;
X    double ehp, ha, dec;    /* ehp: angular dia of earth from body */
X    int new;
X
X    if (same_cir (np, &last_now) && about_now (np, &last_now, as*.00028)) {
X        *sp = last_sky;
X        new = 0;
X    } else {
X        double lsn, rsn;
X        double deps, dpsi;
X
X        last_now = *np;
X        sunpos (mjd, &lsn, &rsn);        /* sun's true ecliptic long
X                         * and dist
X                         */
X        nutation (mjd, &deps, &dpsi);    /* correct for nutation */
X        lsn += dpsi;
X        lsn -= degrad(20.4/3600);        /* and light travel time */
X
X        sp->s_edist = rsn;
X        sp->s_sdist = 0.0;
X        sp->s_elong = 0.0;
X        sp->s_size = raddeg(4.65242e-3/rsn)*3600*2;
X        sp->s_mag = -26.8;
X        sp->s_hlong = lsn-PI;    /* geo- to helio- centric */
X        range (&sp->s_hlong, 2*PI);
X        sp->s_hlat = 0.0;
X
X        ecl_eq (mjd, 0.0, lsn, &last_ra, &last_dec);
X        sp->s_ra = last_ra;
X        sp->s_dec = last_dec;
X        if (epoch != EOD)
X        precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
X        new = 1;
X    }
X
X    now_lst (np, &lst);
X    ha = hrrad(lst) - last_ra;
X    ehp = (2.0 * 6378.0 / 146.0e6) / sp->s_edist;
X    ta_par (ha, last_dec, lat, height, ehp, &ha, &dec);
X    hadec_aa (lat, ha, dec, &alt, &az);
X    refract (pressure, temp, alt, &alt);
X    sp->s_alt = alt;
X    sp->s_az = az;
X    last_sky = *sp;
X    return (new);
X}
X
X/* find moon's circumstances now.
X * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
X * return 0 if only alt/az changes, else 1 if all other stuff updated too.
X */
Xmoon_cir (as, np, sp)
Xdouble as;
XNow *np;
XSky *sp;
X{
X    static Sky last_sky;
X    static Now last_now;
X    static double ehp;
X    static double last_ra, last_dec;    /* unprecessed */
X    double lst, alt, az;
X    double ha, dec;
X    int new;
X
X    if (same_cir (np, &last_now) && about_now (np, &last_now, as*.000021)) {
X        *sp = last_sky;
X        new = 0;
X    } else {
X        double lam, bet;
X        double deps, dpsi;
X        double lsn, rsn;    /* sun long in rads, earth-sun dist in au */
X        double edistau;    /* earth-moon dist, in au */
X        double el;        /* elongation, rads east */
X
X        last_now = *np;
X        moon (mjd, &lam, &bet, &ehp);    /* moon's true ecliptic loc */
X        nutation (mjd, &deps, &dpsi);    /* correct for nutation */
X        lam += dpsi;
X        range (&lam, 2*PI);
X
X        sp->s_edist = 6378.14/sin(ehp);    /* earth-moon dist, want km */
X        sp->s_size = 3600*31.22512*sin(ehp);/* moon angular dia, seconds */
X
X        ecl_eq (mjd, bet, lam, &last_ra, &last_dec);
X        sp->s_ra = last_ra;
X        sp->s_dec = last_dec;
X        if (epoch != EOD)
X        precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
X
X        sunpos (mjd, &lsn, &rsn);
X        range (&lsn, 2*PI);
X        elongation (lam, bet, lsn, &el);
X
X        /* solve triangle of earth, sun, and elongation for moon-sun dist */
X        edistau = sp->s_edist/1.495979e8; /* km -> au */
X        sp->s_sdist =
X        sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
X
X        /* TODO: improve mag; this is based on a flat moon model. */
X        sp->s_mag = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))));
X
X        sp->s_elong = raddeg(el);    /* want degrees */
X        sp->s_phase = fabs(el)/PI*100.0;    /* want non-negative % */
X        sp->s_hlong = sp->s_hlat = 0.0;
X        new = 1;
X    }
X
X    /* show topocentric alt/az by correcting ra/dec for parallax
X     * as well as refraction.
X     */
X    now_lst (np, &lst);
X    ha = hrrad(lst) - last_ra;
X    ta_par (ha, last_dec, lat, height, ehp, &ha, &dec);
X    hadec_aa (lat, ha, dec, &alt, &az);
X    refract (pressure, temp, alt, &alt);
X    sp->s_alt = alt;
X    sp->s_az = az;
X    last_sky = *sp;
X    return (new);
X}
X
X/* given geocentric ecliptic longitude and latitude, lam and bet, of some object
X * and the longitude of the sun, lsn, find the elongation, el. this is the
X * actual angular separation of the object from the sun, not just the difference
X * in the longitude. the sign, however, IS set simply as a test on longitude
X * such that el will be >0 for an evening object <0 for a morning object.
X * to understand the test for el sign, draw a graph with lam going from 0-2*PI
X *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
X *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
X *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
X *   first line and bounded within the second pair of lines.
X * all angles in radians.
X */
Xelongation (lam, bet, lsn, el)
Xdouble lam, bet, lsn;
Xdouble *el;
X{
X    *el = acos(cos(bet)*cos(lam-lsn));
X    if (lam>lsn+PI || lam>lsn-PI && lam<lsn) *el = - *el;
X}
X
X/* return whether the two Nows are for the same observing circumstances. */
Xsame_cir (n1, n2)
Xregister Now *n1, *n2;
X{
X    return (n1->n_lat == n2->n_lat
X        && n1->n_lng == n2->n_lng
X        && n1->n_temp == n2->n_temp
X        && n1->n_pressure == n2->n_pressure
X        && n1->n_height == n2->n_height
X        && n1->n_tz == n2->n_tz
X        && n1->n_epoch == n2->n_epoch);
X}
X
X/* return whether the two Nows are for the same LOCAL day */
Xsame_lday (n1, n2)
XNow *n1, *n2;
X{
X    return (mjd_day(n1->n_mjd - n1->n_tz/24.0) ==
X        mjd_day(n2->n_mjd - n2->n_tz/24.0));
X}
X
X/* return whether the mjd of the two Nows are within dt */
Xstatic
Xabout_now (n1, n2, dt)
XNow *n1, *n2;
Xdouble dt;
X{
X    return (fabs (n1->n_mjd - n2->n_mjd) <= dt/2.0);
X}
X
Xnow_lst (np, lst)
XNow *np;
Xdouble *lst;
X{
X    utc_gst (mjd_day(mjd), mjd_hr(mjd), lst);
X    *lst += radhr(lng);
X    range (lst, 24.0);
X}
X
X/* round a time in days, *t, to the nearest second, IN PLACE. */
Xrnd_second (t)
Xdouble *t;
X{
X    *t = floor(*t*SPD+0.5)/SPD;
X}
X
Xdouble
Xmjd_day(jd)
Xdouble jd;
X{
X    return (floor(jd-0.5)+0.5);
X}
X
Xdouble
Xmjd_hr(jd)
Xdouble jd;
X{
X    return ((jd-mjd_day(jd))*24.0);
X}
EOFxEOF
len=`wc -c < circum.c`
if expr $len != 10181 > /dev/null
then echo Length of circum.c is $len but it should be 10181.
fi
echo x comet.c
sed -e 's/^X//' << 'EOFxEOF' > comet.c
X#include <math.h>
X#include "astro.h"
X
X/* given a modified Julian date, mjd, and a set of heliocentric parabolic
X * orbital elements referred to the epoch of date (mjd):
X *   ep:   epoch of perihelion,
X *   inc:  inclination,
X *   ap:   argument of perihelion (equals the longitude of perihelion minus the
X *       longitude of ascending node)
X *   qp:   perihelion distance,
X *   om:   longitude of ascending node;
X * find:
X *   lpd:  heliocentric longitude,
X *   psi:  heliocentric latitude,
X *   rp:   distance from the sun to the planet,
X *   rho:  distance from the Earth to the planet,
X *   lam:  geocentric ecliptic longitude,
X *   bet:  geocentric ecliptic latitude,
X *         none are corrected for light time, ie, they are the true values for
X *       the given instant.
X *
X * all angles are in radians, all distances in AU.
X * mutual perturbation corrections with other solar system objects are not
X * applied. corrections for nutation and abberation must be made by the caller.
X * The RA and DEC calculated from the fully-corrected ecliptic coordinates are
X * then the apparent geocentric coordinates. Further corrections can be made,
X * if required, for atmospheric refraction and geocentric parallax.
X */
Xcomet (mjd, ep, inc, ap, qp, om, lpd, psi, rp, rho, lam, bet)
Xdouble mjd;
Xdouble ep, inc, ap, qp, om;
Xdouble *lpd, *psi, *rp, *rho, *lam, *bet;
X{
X    double w, s, s2;
X    double l, sl, cl, y;
X    double spsi, cpsi;
X    double rd, lsn, rsn;
X    double lg, re, ll;
X    double cll, sll;
X    double nu;
X
X#define    ERRLMT    0.0001
X        w = ((mjd-ep)*3.649116e-02)/(qp*sqrt(qp));
X        s = w/3;
X    while (1) {
X        double d;
X        s2 = s*s;
X        d = (s2+3)*s-w;
X        if (fabs(d) <= ERRLMT)
X        break;
X        s = ((2*s*s2)+w)/(3*(s2+1));
X    }
X
X        nu = 2*atan(s);
X    *rp = qp*(1+s2);
X    l = nu+ap;
X        sl = sin(l);
X    cl = cos(l);
X    spsi = sl*sin(inc);
X        *psi = asin(spsi);
X    y = sl*cos(inc);
X        *lpd = atan(y/cl)+om;
X    cpsi = cos(*psi);
X        if (cl<0) *lpd += PI;
X    range (lpd, 2*PI);
X        rd = *rp * cpsi;
X    sunpos (mjd, &lsn, &rsn);
X    lg = lsn+PI;
X        re = rsn;
X    ll = *lpd - lg;
X        cll = cos(ll);
X    sll = sin(ll);
X        *rho = sqrt((re * re)+(*rp * *rp)-(2*re*rd*cll));
X        if (rd<re)
X            *lam = atan((-1*rd*sll)/(re-(rd*cll)))+lg+PI;
X    else
X        *lam = atan((re*sll)/(rd-(re*cll)))+*lpd;
X    range (lam, 2*PI);
X        *bet = atan((rd*spsi*sin(*lam-*lpd))/(cpsi*re*sll));
X}
EOFxEOF
len=`wc -c < comet.c`
if expr $len != 2390 > /dev/null
then echo Length of comet.c is $len but it should be 2390.
fi